#include <iostream>
#include <fstream>
using namespace std;

#include "Signal.h"
#include "MPC_EW1D_PatternedW_b.h"

// $project.exe W
const int T_NODE = 1001;
int main(int argc,char** argv){
  double time[T_NODE];
  double value_r2[T_NODE];
  double value_m2[T_NODE];
  double input[2];
  input[0] = atof(argv[1]);
  input[1] = atof(argv[2]);
  ofstream log_r2("log_r2.dat");
  ofstream log_m2("log_m2.dat"); 

  // Setup the model
  double time_to_end  = 0.0;
  double dt           = 1.0;
  int    mode	      = 20000;
  double Xmax	      = 40000;
  int    LatticeSize  = 40000;
  double Freq	      = 0.0;
  
  double* cov_z[2];
  double* mean_z[2];

  cov_z[0]  = new double[mode];
  cov_z[1]  = new double[mode];
  mean_z[0] = new double[mode];
  mean_z[1] = new double[mode];

  for(int m=0;m<mode;m++){
    cov_z[0][m] = 0.0;
    cov_z[1][m] = 0.0;
    mean_z[0][m] = 0.0;
    mean_z[1][m] = 0.0;
  }

  MathFunction *f_r2 = new eq_r2_f(dt,
				   mode,
				   Xmax,
				   LatticeSize,
				   Freq);
  MathFunction *f_m2 = new eq_m2_f(dt,
				   mode,
				   Xmax,
				   LatticeSize,
				   Freq);

  f_r2->parameterize(0,cov_z[0]);
  f_r2->parameterize(1,mean_z[0]);
  f_r2->parameterize(2,cov_z[1]);
  f_r2->parameterize(3,mean_z[1]);

  f_m2->parameterize(0,cov_z[0]);
  f_m2->parameterize(1,mean_z[0]);
  f_m2->parameterize(2,cov_z[1]);
  f_m2->parameterize(3,mean_z[1]);

  for(int i=0;i<T_NODE;i++){
    time[i]=double(i);

    f_r2->parameterize(4,time[i]);
    value_r2[i] = f_r2->value(input);
    log_r2 << setw(8) << time[i] << setw(15) << value_r2[i] << setw(3) << 0.0; // The last 0.0 is only to make the output has 3 columns
    if(i!=T_NODE-1){
      log_r2 << "\n";
    }

    f_m2->parameterize(4,time[i]);
    value_m2[i] = f_m2->value(input);
    log_m2 << setw(8) << time[i] << setw(15) << value_m2[i] << setw(3) << 0.0;
    if(i!=T_NODE-1){
      log_m2 << "\n";
    }
  }
  delete f_r2;
  delete f_m2;
  delete cov_z[0];
  delete cov_z[1];
  delete mean_z[0];
  delete mean_z[1];
}